function [X_cont_ix_tod,X_cont_ix,tod] = fn_cont_ix_tod(tr_alpha,tr_beta,Xret,deltan,nDays,nObs_per_day,tod)
% IN
%      tr_alpha         jump truncation para
%      tr_beta
%      Xret             size [p, n], and this allows for p>1
%      deltan
%      nDays            number of DAYS, cannot be other time intervals.
%      nObs_per_day
% OUT
%      X_cont_ix_tod    logicals: is Xret cont? size [p n]. This version
%      adjusts the diurnal effect.
%      tod              Estimates of the diurnal effect.
%      X_cont_ix        logicals: is Xret cont? size [p n]. No adjustment
%      for diurnal effect.

p      = size(Xret,1);

if isempty(tod)
    avgvol = zeros(p,1);
    tod    = ones(p, nObs_per_day);


    if isvector(Xret)
 
        X_cont_ix = fn_jump_ix_1variable_tod(tr_alpha,tr_beta,Xret,deltan,nDays,nObs_per_day,tod);

    else
  
        X_cont_ix = -999*ones(size(Xret));

        for idx_p = 1:p
    
            X_cont_ix(idx_p,:) = fn_jump_ix_1variable_tod(tr_alpha,tr_beta,...
                                          Xret(idx_p,:),deltan,nDays,nObs_per_day,tod(idx_p,:));
        end

    end

    for j = 1:p
    
        avgvol(j) = sum((Xret(j,:).^2).*X_cont_ix(j,:));

        for i = 1:nObs_per_day

            tod(j,i) = nObs_per_day*sum((Xret(j,i:nObs_per_day:end).^2).*(X_cont_ix(j,i:nObs_per_day:end)==1));

        end

    end

    tod = (tod./repmat(avgvol,1,nObs_per_day));

end

if isvector(Xret)
 
    X_cont_ix_tod = fn_jump_ix_1variable_tod(tr_alpha,tr_beta,Xret,deltan,nDays,nObs_per_day,tod);

else
  
    X_cont_ix_tod = -999*ones(size(Xret));

    for idx_p = 1:p
    
        X_cont_ix_tod(idx_p,:) = fn_jump_ix_1variable_tod(tr_alpha,tr_beta,...
                                          Xret(idx_p,:),deltan,nDays,nObs_per_day,tod(idx_p,:));
    end

end





function Xix = fn_jump_ix_1variable_tod(tr_alpha,tr_beta,Xret,deltan,nDays,nobs,tod)
trunc_const = tr_alpha.*deltan^(tr_beta);        
XretMat     = reshape(Xret,nobs,nDays);
var1         =   (sum(XretMat(1:end,:).^2)*252);
var2         =   sum(XretMat(1:end,:).^2.*(abs(XretMat(1:end,:))<trunc_const*sqrt(var1)))*252;
var2         =   kron(var2,tod);
Xix          =   abs(Xret)<trunc_const*sqrt(var2);